var airpv = [];
var airvx = [];
var airvy = [];
var outpv = [];
var outvx = [];
var outvy = [];
var blockair = [];

for (var i = 0; i < 50; i++)
{
	airpv [i] = new Float64Array(50);
	airvx [i] = new Float64Array(50);
	airvy [i] = new Float64Array(50);
	outpv [i] = new Float64Array(50);
	outvx [i] = new Float64Array(50);
	outvy [i] = new Float64Array(50);
	blockair [i] = new Int8Array(50);
}

var kernel = [];

(function ()
{
	var i, j;
	var s = 0.0;
	for (j=-1; j<2; j++)
		for (i=-1; i<2; i++)
		{
			kernel[(i+1)+3*(j+1)] = Math.exp(-2.0*(i*i+j*j));
			s += kernel[(i+1)+3*(j+1)];
		}
	s = 1.0 / s;
	for (j=-1; j<2; j++)
		for (i=-1; i<2; i++)
			kernel[(i+1)+3*(j+1)] *= s;
})();

function air_sim ()
{
	// define AIR_PLOSS 0.9999
	// define AIR_TSTEPP 0.3
	// define AIR_VLOSS 0.999
	// define AIR_TSTEPV 0.4
	// define AIR_VADV 0.3

	var i, j, ii, jj, x, y;
	var dp, dx, dy, f, tx, ty;
	var advDistanceMult = 0.7;
	var stepX, stepY;
	var stepLimit, step;

	for (i = 0; i < 50; i++) // Pressure is decreased at the edges
	{
		airpv[i][ 0] *= 0.8;
		airpv[i][ 1] *= 0.8;
		airpv[i][ 2] *= 0.8;
		airpv[i][48] *= 0.8;
		airpv[i][49] *= 0.8;

		airvx[i][ 0] *= 0.9;
		airvx[i][ 1] *= 0.9;
		airvx[i][48] *= 0.9;
		airvx[i][49] *= 0.9;

		airvy[i][ 0] *= 0.9;
		airvy[i][ 1] *= 0.9;
		airvy[i][48] *= 0.9;
		airvy[i][49] *= 0.9;
	}

	for (i = 0; i < 50; i++) // Pressure is decreased at the edges
	{
		airpv[ 0][i] *= 0.8;
		airpv[ 1][i] *= 0.8;
		airpv[ 2][i] *= 0.8;
		airpv[48][i] *= 0.8;
		airpv[49][i] *= 0.8;

		airvx[ 0][i] *= 0.9;
		airvx[ 1][i] *= 0.9;
		airvx[48][i] *= 0.9;
		airvx[49][i] *= 0.9;

		airvy[ 0][i] *= 0.9;
		airvy[ 1][i] *= 0.9;
		airvy[48][i] *= 0.9;
		airvy[49][i] *= 0.9;
	}

	for (j = 1; j < 50; j++) //clear some velocities near walls
	{
		for (i = 1; i < 50; i++)
		{
			if (blockair[j][i])
			{
				airvx[j][i] = 0;
				airvx[j][i-1] = 0;
				airvy[j][i] = 0;
				airvy[j-1][i] = 0;
			}
		}
	}
	
	for (y = 1; y < 50; y++) //pressure adjustments from velocity
	{
		for (x = 1; x < 50; x++)
		{
			dp  = airvx[y][x-1] - airvx[y][x];
			dp += airvy[y-1][x] - airvy[y][x];
			airpv[y][x] *= 0.9999; // AIR_PLOSS
			airpv[y][x] += dp * 0.3; // AIR_TSTEPP
		}
	}
	
	for (y = 0; y < 49; y++) //velocity adjustments from pressure
	{
		for (x = 0; x < 49; x++)
		{
			dx = airpv[y][x] - airpv[y][x+1];
			dy = airpv[y][x] - airpv[y+1][x];
			airvx[y][x] *= 0.999; // AIR_VLOSS
			airvy[y][x] *= 0.999; // AIR_VLOSS
			airvx[y][x] += dx * 0.4; // AIR_TSTEPV
			airvy[y][x] += dy * 0.4; // AIR_TSTEPV
			if (blockair[y][x] || blockair[y][x+1])
				airvx[y][x] = 0.0;
			if (blockair[y][x] || blockair[y+1][x])
				airvy[y][x] = 0.0;
		}
	}
	
	for (y=0; y<50; y++) // update velocity and pressure
	{
		for (x=0; x<50; x++)
		{
			// The air sim code does use only the previous frame to generate the next frame
			dx = 0.0;
			dy = 0.0;
			dp = 0.0;
			for (j=-1; j<2; j++)
			{
				for (i=-1; i<2; i++)
				{
					f = kernel[i+1+(j+1)*3];
					jj = y+j; ii = x+i;
					if (jj < 1 || jj > 48 || ii < 1 || ii > 48 || blockair[jj][ii])
					{
						jj = y; ii = x;
					}
					dx += airvx[jj][ii] * f;
					dy += airvy[jj][ii] * f;
					dp += airpv[jj][ii] * f;
				}
			}
			tx = x - dx * 0.7;
			ty = y - dy * 0.7;
			
			if ((Math.abs(dx * 0.7) > 1 || Math.abs(dy * 0.7) > 1) && (tx >= 2 && tx < 48 && ty >= 2 && ty < 48))
			{
				// Trying to take velocity from far away, check whether there is an intervening wall. Step from current position to desired source location, looking for walls, with either the x or y step size being 1 cell
				if (Math.abs(dx) > Math.abs(dy))
				{
					stepX = (dx < 0.0) ? 1 : -1;
					stepY = -dy / Math.abs(dx);
					stepLimit = (Math.abs(dx * 0.7)) | 0;
				}
				else
				{
					stepY = (dy < 0.0) ? 1 : -1;
					stepX = -dx / Math.abs(dy);
					stepLimit = (Math.abs(dy * 0.7)) | 0;
				}
				tx = x;
				ty = y;
				for (step = 0; step < stepLimit; ++step)
				{
					tx += stepX;
					ty += stepY;
					if (blockair[(ty+0.5)|0][(tx+0.5)|0])
					{
						tx -= stepX;
						ty -= stepY;
						break;
					}
				}
				if (step === stepLimit)
				{
					// No wall found
					tx = x - dx * 0.7;
					ty = y - dy * 0.7;
				}
			}
			i = tx | 0;
			j = ty | 0;
			tx -= i;
			ty -= j;

			if (!blockair[y][x] && i>=2 && i<=47 && j>=2 && j<=47)
			{
				dx *= 1.0 - 0.3;
				dy *= 1.0 - 0.3;

				dx += 0.3*(1.0-tx)*(1.0-ty)*airvx[j][i];
				dy += 0.3*(1.0-tx)*(1.0-ty)*airvy[j][i];

				dx += 0.3*tx*(1.0-ty)*airvx[j][i+1];
				dy += 0.3*tx*(1.0-ty)*airvy[j][i+1];

				dx += 0.3*(1.0-tx)*ty*airvx[j+1][i];
				dy += 0.3*(1.0-tx)*ty*airvy[j+1][i];

				dx += 0.3*tx*ty*airvx[j+1][i+1];
				dy += 0.3*tx*ty*airvy[j+1][i+1];
			}
			
			outvx[y][x] = dx;
			outvy[y][x] = dy;
			outpv[y][x] = dp;
		}
	}
	
	for (y=0; y<50; y++) // update velocity and pressure
	{
		for (x=0; x<50; x++)
		{
			airvx[y][x] = outvx[y][x];
			airvy[y][x] = outvy[y][x];
			airpv[y][x] = outpv[y][x];
		}
	}
}